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Functioning as sensors and propulsors, cilia are evolutionarily conserved organelles having a highly 
organized internal structure. How a Paramecium's cilium produces off-propulsion-plane curvature during 
its return stroke for symmetry breaking and drag reduction is not known. We explain these cilium 
deformations by developing a torsional pendulum model of beat frequency dependence on viscosity and an 
olivo-cerebellar model of self-regulation of posture control. The phase dependence of cilia torsion is 
determined, and a bio-physical model of hardness control with predictive features is offered. Crossbridge 
links between the central microtubule pair harden the cilium during the power stroke; this stroke's end is a 
critical phase during which ATP molecules soften the crossbridge-microtubule attachment at the cilium 
inflection point where torsion is at its maximum. A precipitous reduction in hardness ensues, signaling the 
start of ATP hydrolysis that re-hardens the cilium. The cilium attractor basin could be used as reference for 
perturbation sensing. 



The paramecium is slipper-shaped ciliate protozoa widely found in oxygenated aquatic environments. These 
animals propel themselves, albeit with limited maneuverability, by the synchronous motion of numerous 
tiny cilia populated around their flexible bodies. Paramecia are 100 - 350 urn long, deformable, and may 
contain up to 3000 flexible cilia, each being about 17 um long and 0.25 um in diameter 1 . For propulsion, each 
cilium beats at about 14.1 Hz in water while undergoing a complicated, three-dimensional, phase-dependent 
motion. Cilia have a highly organized internal structure 2 . The assemblage of cilia also undergoes a coordinated 
motion, called metachronic motion — i.e., the assemblage beats with a constant phase difference between neigh- 
boring cilia 3 . High-speed stereomicroscopic measurements of cilia motion show that during the power stroke, the 
entire cilium — from the root to the distal point — remains virtually in a plane in a straight and taut posture 4 . The 
return stroke, on the other hand, is a continuously deforming, three-dimensional motion. Near the end of the 
power stroke, the cilium abruptly wilts to a low-drag posture, initiating the return stroke with a shift to the lateral 
plane. During the return stroke, which initially lies in a plane parallel to the base surface, the cilium gradually 
becomes straight, erecting to an upright posture while approaching the start of the power stroke. In the tiny 
Paramecium, the posture of each of its 3000 cilia switches between taut and slack every time they beat — in 0.07 s — 
a remarkable feat. Considering the metachronism property, the temporal accuracy and robustness of the phase 
control of posture between the neighboring cilia is also extraordinary. If the change in posture is due to hardness 
control, then how does the cilium's hardness abruptly switch between high and low values? (Hardness is a 
material property given by the moduli of elastic and torsional rigidity.) This switching in posture has not been 
explained by any theory. Neither is there any reported direct experimental or biochemical evidence of any 
hardness control. Here, we consider theoretical modeling of the cilium deformation, proposing a new mechanism 
for hardness control. 

Understanding of the mechanism for controlling hardness may help us build novel propulsive actuators that 
interact with vortex flows to absorb their energy and release it at appropriate phase 5-6 . Some of these cilia act as 
sensors 7,8 , and an understanding of hardness control may help us build (currently unavailable) nonlinear sensors. 
In the cilia literature, we have not found any awareness of hardness control; in current nanorod fabrications of 
biomimetic cilia, the hardness is not controlled 9 ; the cilium motion is approximated to be two-dimensional in 
notable theories on metachronism (torsion is ignored) 10 . In the context of research on the emerging multisystemic 
disorders called ciliopathy, the present work defines one facet of the interaction between the structure-function 
relationship and signaling pathways in a normal cilium 7,11 . 
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Below, we present four models. Overall, assuming energy conser- 
vation to be of paramount importance, in the early models we estab- 
lish that drag is determined by a self- regulating process; resonant 
amplified oscillations (low friction) take place; and the cilium has 
locked-in (same initial condition) nonlinear (repeatable) trajectories; 
then the resulting phase relationship is used to develop the final 
breakup-makeup model. In the first model, the balance between 
the applied moment at the cilium base and fluid drag on the cilium 
is considered. The cilium base, which is partly constrained, experi- 
ences a circumferential perturbation actuated by a helicopter swash 
plate-like mechanism 2 . The distal end, however, is free, and there is a 
time lag in response between the two ends. As a result, the cilium is 
subjected to torsion and is modeled as a torsional pendulum. In the 
second model, the pendulum is shown to have a nonlinear character. 
The critical phase when the pendulum is at neutral equilibrium (the 
velocity-phase map is a saddle) is identified and its significance is 
discussed. In the first part of the third model, the position data of the 
cilium from the literature are digitized and processed to unambigu- 
ously identify the beginning and end of the power stroke. In the 
second part of this model, spatio-temporal curvature and torsion 
maps of the cilium are determined. This exercise helped to identify 
the phase of the beat cycle and the location on the cilium where 
torsion reaches its maximum value — that is, where the most vulner- 
able point for structural failure lies. In the fourth model, the mech- 
anism of how the hardness is switched at this point to result in a 
transition in the cilium from taut to wilting to taut again is given. 
Efficiency and the precision of damping control are discussed. 
Finally, the theoretical foundation of nonlinear sensing, in the con- 
text of the formulation of self-regulation, is discussed. 

Results 

Beat frequency is inversely proportional to the square root of vis- 
cosity. First, we carry out a simple modeling of the measurements of 
how the beat frequency of the Paramecium's cilium varies with the 
absolute viscosity p., which varies as the reciprocal of water tem- 
perature 3 . This variation is shown in Fig. 1. The cilium's Reynolds 



number (Re) — the ratio of inertia to viscous forces — is extremely low 
(Re= (25-875) x 10~ 6 ) if the 100- to 350-um-long Paramecium 
moves 1 to 10 body lengths per second, the cilia are 0.25 urn in 
diameter (d), and kinematic viscosity of water (u) is 1 x 10 ~ 6 . This 
low Reynolds number range is called the Stokes flow regime. We 
assume here that Re is varied due to variations in v (= pip, where 
p is the density of water) only and that the velocity and length scales 
remain constant. In other words, the Paramecium is assumed to 
always move at the same speed, forward or backward. (While the 
Paramecium has a limited portfolio of maneuverability, it is shown 
later that it has full self- regulation over the motion it is capable of.) 
We assume that the cilium has vanishingly small frictional resistance 
at its base and is undergoing resonant amplified oscillation. In other 
words, the quality factor (Q) is very high 12 . The cilium is just under 
critically damped (damping ratio f < 1.0) so that it is ringing but 
would stop if raised to £ = 1.0 — in other words, it is completing the 
power and return strokes in almost minimum time. There is no 
undue oscillation at the end of the strokes. These assumptions 
imply that energy is optimally consumed; therefore, there are 
robust ionic mechanisms (capable of disturbance rejection) to 
ensure that the damping is kept a bare tad below what is required 
for critical damping. 

We assume that the cilium acts like a torsional pendulum. Energy 
is spent to wind it like a torsional spring during the return stroke, and 
the energy is released during the natural unwinding during the power 
stroke, thus producing thrust. Conversely, if the cilium acts like a 
sensor, then the extracellular perturbation produces ionic pulses at 
the cilium's base. As shown in Fig. la, at this time we disregard the 
shape of the cilium and assume that a torque x is applied at the base 
and the cilium experiences a drag torque (D T ) in the opposite 
direction. 

For small rotational angles 0, x = — z k0, where t k is the torsional 
spring constant. If / is the moment of inertia (which remains con- 
stant), the simplified equation of motion is JO = x = T k 0. If/ D (1/s) is 
the cilium beat frequency, the angular beat frequency is a) 0 = 2nf„, 
and the solution is co 0 = \J z kj]. We assume that T kccRecc l/p, that 
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Figure 1 | Torsional modeling of the variation of cilium beat frequency with viscosity, (a) Schematic of the torsional pendulum model of the cilium; 
(b) Comparison of the model (solid line: f 0 = K(l/J(i), K = 75, and ji in cp) with measurements 3 (dots). 



SCIENTIFIC REPORTS | 3 : 1956 | DOI: 10.1038/srep01956 



2 



is, the torsional sp ring constant is inversely proportional to viscosity. 
Therefore, f 0 oc yj\/f.i. This solution is compared with measure- 
ments 3 in Fig. lb a nd th e agreement is good. (Conversely, since the 
data follows/,, cc a/1/ H, T kccReaz 

This model gives us the hypotheses that the cilium may be treated 
as a resonant torsional pendulum. We next model the trajectory 
of the cilium pendulum and explore if it behaves like a nonlinear 
pendulum. 

Cilium motion follows olivo-cerebellar dynamical systems 
equations. In this section, olivo-cerebellar dynamics, which owes 
its origin to inferior-olive (10) neuron dynamics 1319 , is applied to 
the cilium motion. Voltage stimulation affects cilia motility 3,20,21 . 
Slow voltage ramps (5 to 10 mV/s) and moderate amplitudes 
(±20 mV) make the membrane properties time dependent, which 
affects ciliary response (cycle time and relaxation time). The cilium 
beat frequency and reverse/forward motion switching are affected by 
calcium channel activation 1,22 . The alignment of ciliary beating to the 
fluid flow direction is controlled by a mechanism in which certain 
cilia act as sensors for lowpass filtering of hydrodynamic forces and 
for generating a polarity signal for directional alignment 23 . There is 
direct evidence of cilium beat frequency control by Ca 2+ -ions 21 . An 
examination of 60 data sets on flagella and cilia shows that stretch- 
activated calcium channels cause calcium and other waves to 
propagate at speeds of 100 to 1000 um/s, which affects the tubular 
actuators 24 . We assume that Ca 2+ -ions control cilium motion 
(describable as a primitive neuron 14 ). It is also noted that IO 
neurons are merely one of many situations in nature where Ca 2+ - 
oscillator type equations such as that below apply 5,6 . 
The model of the i th ion-related controller 171 " is given as 
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where the variables z,- and w, are associated with sub-threshold oscil- 
lations and low-threshold (Ca 2+ -dependent) spiking in olivo- 
cerebellar dynamics (the higher threshold (Na + -dependent) spiking 
oscillator is not considered). Here, equation (1) is referred to as the 
Ca 2+ -oscillator equation (FitzHugh-Nagumo model). The constant 
parameter e, Ca controls the oscillation time scale; and I Ca drives the 
depolarization levels. The nonlinear function is p, z (z,-) = Z;(z,— a,-) 
(1 — Zj), where a,- is a constant parameter associated with the non- 
linear function. The function I ext (f) is the extracellular stimulus, 
whose amplitude and duration are used for the purpose of control 
(changing the motion of Paramecium direction, for example). If I ext 
(r) = 0 (independent oscillator), equation (1) can be written as 
Zi + F(zi)zi + kz, +el = 0, where F is a cubic polynomial function 
and A: is a constant. This equation bears resemblance to Lienard's 
oscillator (in contrast, the function F is a well-defined quadratic in 
the familiar van der Pol oscillator) 25,26 . The oscillator exhibits a 
closed-orbit T, in the state space (z, — Z;), that is, (z, — w,-), which 
is also known as limit cycle oscillation (LCO), the constant para- 
meters determining the form of J 1 ,-. 

The Cfl 2+ -oscillator equation is solved using an electronic analog 
oscillator 27 . The (z, w) waveform data are gathered from Simulink 
and saved into files for reading; these are z, w waveforms, their 
derivatives, and their second derivatives (see Methods and Supple- 
mentary Information (SI. 1)). In Fig. 2, the phase maps are scaled to 
the same size as the experimental data and centered at the same point 
so that the waveforms can be compared. The modeled phase maps 
are also rotated to the proper orientation, which is a 10° rotation for 
position data (z, w) and —55° for derivative data (z,w). (These are 
linear translations of unit V to unit um, and hence should be per- 
missible. This is similar to the calibration of the placement of inertial 
navigation units with respect to the center of gravity (or pressure) in a 
manmade swimming or flying machine.) 



0.3 r 1 




Figure 2 | Comparison of olivo-cerebellar dynamics model with presently 
processed cilium measurements, (a) Position (um/s) (z, w), and (b) 
velocity (um/s 2 ) (z,w). Blue: model; red: presently processed 
measurements of cilium distal point. Parameter a;. 0.015. The numbers 
denote the time steps (TS) in the cilium (see Methods). 



The modeled states (z, w) are compared with measurements in 
Fig. 2 (the trajectories are self-similar along the cilium length (Fig. 6 
in Ref. 3)). Note that the experimental data are for axes shown in 
Figs. 6d and 6e for the cilium for a given value of s. The agreement, 
even in the first derivative, is good for a, = 0.015. The simulation 
shows that the cilium follows a track that is describable as a limit 
cycle. Therefore, it has an autonomic character; the cilium is con- 
trolled without sensors. The degree of control exhibited by the self- 
similarity of the cilium position and velocity all along its length, and 
through the stroke transitions, is remarkable; the present model 
relates the postural self-similarity to a self-regulating process of ionic 
control especially through the stroke transitions. 

The (z = 0,w = 0)location is an inflection point (Fig. 2b) in both 
the biological cilium and the modeled state map (as a result, w is very 
high at TS14; at this phase, as shown in Figs. 3f and 4b, torsion 
reaches a local maximum value at the point marked P on the cilium, 
which lies approximately half way along its length). The significance 
from the point of view of the nonlinear oscillatory behavior of the 
cilium is that point P occurs during a phase where the cilium is at 
neutral equilibrium and the phase map is a saddle. From the design 
point of view, small perturbations applied at this phase will amplify 
rapidly — an ideal strategy for minimizing energy consumption. The 
implication is that this is the best phase to transition from the power 
to the return stroke to maximize thrust and minimize drag, allowing 
the transition to be rapid. If this is accepted, then it only remains for 
us to determine the exact location along the length of the cilium 
where its hardness switches from hard to soft. Surely, it is not prudent 
to spend energy to soften the entire cilium if we could wilt it by 
softening it at one critical location. Therefore, having found the 
critical phase in the beat cycle, we now search for this critical location 
along the length of the cilium and determine the mechanism of soft- 
ening and hardening of the cilium. For a Paramecium, which has to 
execute this hardness transition in so many cilia and so many times 
every second over its entire life, this mechanism must be very low 
cost. 

This model gives us two hypotheses: the cilium switches from hard 
to soft, and that the switching occurs at the saddle point. 

When does the power stroke end? The end of the power stroke (and 
its beginning) can be ambiguous because the proximal part — where 
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Figure 3 | Properties of deformation in the cilium. (a) Local axial angle of the orbital velocity, and (b) local angle of the position of the cilium length 
segments with respect to the x-axis (power stroke axis) during the beat cycle (the graph is independent of cilium length (50 or 17 urn)). The horizontal 
lines in (a) mark angle values of 45° and 90°. (c) Local cilium velocities (um/s); the numbers denote time steps, (d) Local cilium acceleration (blue arrows: 
um/s 2 ). Color code in a to c: different color signifies every 5% of cilium length starting from the base, (e) Components of acceleration (um/s 2 ). 
(f) Average cilium radius (R av ^) of the cilium during the beat cycle (P is a point of inflection on the cilium at TS14; see Fig. 4b.) (g) Equivalent force (or 
power) variable F D (Vd eq measured by um 2 /s) in the cilium, indicating that the stroke boundaries lie where the product crosses the zero value. 




Figure 4 | Curvature and torsion in the cilium. (a) Curvature, (b) torsion. P is a point of inflection on the cilium (see Fig. 5).) 
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the moment is applied — starts the return stroke prior to the distal 
part, which is free. To accurately determine the extent of the power 
stroke, past cilium data were carefully digitized and processed (see 
Methods and SI.2). The properties of the 17-umTong cilium are 
shown in Fig. 3. Figs. 3a and 3b show the variation of position 
and velocity angles subtended to the propulsion direction at 5% 
steps along the length of the cilium (20 values of x ve i ocity ). 

These angles are given by (acos^ len & h ^j * 180/rc^ for position, 



and 



/ X V elocity\ 180 

acos * 

V Vtotal J K 



for velocity. Here, Xi ength , is the 



component of cilium position in three-dimensional space, x ve i ocity 
is the x-component of velocity, and V tota i is the total magnitude of 
all (x, y, z) components of velocity. There are two regions of 
organized behavior, one where the angles are invariant and 
another where they vary. Figs. 3c and 3d show local velocities and 
accelerations, respectively, in three dimensions. The properties are 
organized. The velocity field folds at (x = 0,_y = 0,z = 0) (a point of 
inflection) (z~0 between TS1 and TS16). This time instant occurs 
nearly at the shift from the power stroke to the return stroke. As per 
the nonlinear modeling given earlier, this is the phase of neutral 
equilibrium in a nonlinear pendulum (i.e., the vertical position in a 
nonlinear pendulum if the string is replaced by a solid rod; the 
topology is a saddle — a neutral equilibrium position that is very 
sensitive to disturbances; 180° later, it is a stable focus). The entire 
cilium has an inward acceleration. Fig. 3d shows that the cilium 
induces an inward "tornado," with the magnitude increasing 
distally from the base. (Beyond the distal point, viscous (thrust) jet 
streams are produced 9 .) The literature has several two-dimensional 
hydrodynamic theories, but two-dimensional cilium motion cannot 
produce the tornado. The acceleration limits drop by a factor of 3 if 
the cilium length is reduced from 50 to 17 urn. To compare, in a rigid 
biorobotic cilium, tip speeds of 700-800 um/s and tracer speeds of 
200 um/s have been reported 9 . 

The components of acceleration are shown in Fig. 3e. Past 
research 4 does not explicitly define the extent of the power stroke, 
but it suggests that the stroke extends from where x- acceleration 
crosses the zero value to where it reaches the maximum value. 
Defining the power stroke as when the cilium is moving in the 
negative propulsion direction seems most reasonable, thereby push- 
ing the cilium in the positive propulsion direction. This is also when 
acceleration is increasing in Fig. 3e — approximately from TS35 to 
TS14. The total duration of the power stroke is approximately 0.50 of 
the beat period, which is longer than the reported range of 0.14 to 
0.28 4 . 

The extent of the power stroke is modeled in another manner. An 
average radius of the cilium (R avg ) was defined as follows. The aver- 
age cilium radius (curved length projected on the transverse plane 

y-z) at a given time instant is i? B tg(f) = \J {^z(0} 2 /2» where l yz (t) is 

the curved length of the entire cilium in the y-z plane (transverse 
plane) (sum of the y- and z-components of each length segment of 
the cilium) at a given time instant. To clarify, R avg is not a radial 
distance from the origin to the point on the cilium where R avg lies. 
(The pectoral fins of large swimming and flying animals similarly 
oscillate in a plane transverse to the direction of propulsion 6 .) The 
motion of the cilium is modeled as a (helicopter) disc whose velocity 
varies radially 11 . The disc area is divided equally at R avg such that the 
momentum imparted to the fluid by the inner and the outer areas is 
equal. The variation of R avg with phase is shown in Fig. 3f. The R avg 
values (these are curved length radius averages) were multiplied by 
the local velocity in the propulsion direction (x-axis) to determine 
the directionality and relative scale of the momentum imparted to 
the fluid. At TS15, R avg is a minimum; P is a point of inflection on the 
cilium at TS14, just prior to when R avg is a minimum. As defined 
below, the equivalent force variable (or power) F D {FoccV avg R avg ) 



crosses zero at TS35 and TS14 (Fig. 3g). This is an indication that the 
power stroke extends from TS35 to TS14, and the return stroke 
extends from TS14 to TS35. TS14 is a boundary in the control law, 
x- acceleration, F D and torsion K (shown later in Fig. 4b). Figs. 3f and 
3g show that the reduction of R avg (the reduction in the hardness 
causing a transverse "unfurling" of the cilium, as discussed later) is a 
mechanism for reducing drag during the cilium's return stroke; con- 
versely, the full extension of R avg (increase in hardness) is a mech- 
anism for maximizing output power during the power stroke. The 
x- acceleration in Fig. 3e has maxima and minima at TS14 and TS35, 
respectively — the boundaries of the strokes. The phase where the 
power stroke ends and the return stroke resumes plays a key role 
in the breakup-makeup model given later. 

The use of R avg has produced consistent results and is used below 
to estimate hydrodynamic (thrust) efficiency. Drag force, assumed to 

1 



apply at R t 



av P is given by F D = Q - p V avg 



A, where V avx is the cilium 



velocity at R avg , Cd is the drag coefficient, and d is the diameter of the 
cilium. The projected area of the cilium is A = 2R avg d, where by 

FD = Cd-pV avg 2 2R avg d. If for Stokes flow, Cd = C* /Red, where C* 

is a constant and Re d is the Reynolds number based on cilium dia- 
meter {Red — V avg d/v), then F£}OzV avg R avs . Applying a mean ad- 
vance to the base (= the velocity of the Paramecium, which 
cancels out below, whereby F D is equivalent to power), the hydro- 

dynamic efficiency is given by rj = '- — 



\VR avg \dt 



Jl\VR avg \dt 

where PS and RS are the power and return strokes, respectively, 
and T is the time period of cilium oscillation. Considering the area 
about the zero crossings, the thrust efficiency (the difference between 
the positive and negative areas divided by the total area under the 
curve) is estimated to be 0. 125. This estimate is used in the Discussion 
section to show that the oxygen budget for thrust is very small, thus 
indicating optimization. 

Cilium wilts where torsion peaks ending power stroke. To 

understand the properties of three-dimensional deformation, first 
a two-dimensional hydrodynamic modeling of the cilium motion 
was carried out (see section SI. 3). The cilium was divided into 
numerous links connected by ball joints, where bending rigidity 
was lumped. Cilium stiffness was increased during the power 
stroke and relaxed during the return stroke, and a mean speed-of- 
advance was added to the base motion. Both the rotational motion at 
the base and the drag force at the joints were included. Even in two 
dimensions, the asymmetric beating pattern (Fig. ST4) and the mean 
advance of the base allow for a non-zero thrust (Fig. ST5) obtained 
from the integral of the base force over a beat cycle. The accuracy of 
the two-dimensional model is unknown due to the non-availability 
of thrust and rigidity measurements for the biological cilium. For the 
same reasons, an extension to three dimensions did not seem very 
productive. However, the three-dimensional model of cilium motion 
allows out-of-plane curvature and torsion, but the two-dimensional 
model does not. 

Consider the curvature fc and torsion t of the cilium. Define r (s, r) 
as the position vector of a point on the surface of the cilium. Primes 
denote the derivatives with respect to s, x denotes the cross product, 
and 1 1 is magnitude. The following equations were used to calculate k 

r' x v"\ 

and x for the cilium data 28 . Curvature is defined as K = 5 — , and 



torsion is defined as t = 



\r'xr"\ 2 



k measures the deviance of a 



curve from being a straight line relative to the osculating plane, and z 
measures how sharply it is twisting. For x = 0, the curve lies com- 
pletely in the same osculating plane (there is only one osculating 
plane). Note that K and t of a helix are constant; when they are not 
constant, the geometry is not helical as in a cilium as opposed to a 
spermatozoon. 
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(b) 

Figure 5 | Schematic representation of the breakup and then makeup mechanics during the transition from the power to the return strokes. 

(a) Schematic of compression-tension neutral-plane flipping at point P on the cilium at the beat phase of neutral equilibrium (see inset P). Symbols: 9, 14, 
and 16 are time steps; m and -m are moments in the x-y plane at the base; D, V, and r are drag, velocity, and the radius of curvature; the prime sign (') 
denotes values near the base; C and T are compressive and tensile stresses; and N = neutral plane. The cilium is shifted along X with each time step to avoid 
clutter. Left inset: CM al4 is the central microtubule a at TS14; CM bl4 is the other central microtubule of the pair at TS14; CB P14 is the crossbridge link 
(CBL) at the point of inflection P at TS14. The CBL can be expected to fail first at the point marked by the red dot in inset P. (b) Modeling of the CBL 
between the pair of central microtubules (the 2 of the 9 + 2). Only the CBL at point P is shown; this is the only one that cracks because torsion reaches high 
values at this point; other CBLs add to the hardness, but they do not crack because the torsion is not high (the region near the cilium tip is ignored) . Broken 
lines represent later times. Curvature is small near point P at TS14; therefore, planar vibrations have a negligible effect on cracking. 



The values of K and z in the cilium are compared in Fig. 4. The 
white lines indicate the boundaries of the power and return strokes; 
the return stroke lies between the white lines. TS1 and TS42 of the 
cilium are assigned a phase of 0 and 27t, respectively; that is, the phase 
at the nth time step is (j> = 2nn/N, where N (= 41) is the total number 
of time steps. The results are lowpass filtered. 

The simulations confirm the past cilium result 4 that k is small 
during the power stroke (because the cilium is stiff), and that K 
and i return to minimal values by the beginning of the power stroke. 
(This finding can be further verified by looking at the position data in 
Fig. 6f.) In the cilium, curvature is high near the base after the return 
stroke has started. Torsion reaches a peak approximately halfway 
along the length when the power stroke ends (marked as point P 
in Fig. 4b). The limits of k and t widen if the cilium length is reduced 
from 50 to 17 p:m. 

This section shows that torsion reaches a peak value halfway along 
the length of the cilium at the saddle point of the stroke phase, which 
was identified earlier using FitzHugh-Nagumo control model. 

Breakup and makeup model of cilium hardness control. The cur- 
vature and torsion results are synthesized in Fig. 5 to understand the 
behavior at the transition from the power to the return stroke when 
the hardness of the cilium changes precipitously. 

The model schematic in Fig. 5a shows projections of the cilium 
onto the propulsion plane at TS9, TS14, andTS16. (These projections 
are copies from Fig. 6f.) Between TS13 and TS14, the time derivative 
of the angle subtended by the cilium at the base changes sign (the 



cilium positions converge in Fig. 6f). This inference is depicted in 
Fig. 5a by the change in the sign of the moments applied at the base 
m. The directions of velocities (V) of the cilium and of the drags (D) 
acting on it are shown. At TS14, while the lower part of the cilium 
turns counterclockwise, the distal part turns in the opposite dir- 
ection. As a result, a point of inflection is produced in between 
(marked as P), where the finite radius of curvature (r) changes sign. 
The cilium precipitously droops after TS14, which indicates a 
decrease in hardness. The cilium is straight during the power stroke. 
Therefore, the question is: How does the hardness drop and recover 
in every beat cycle? This process is modeled in insets P, PI, and P2 in 
Fig. 5a and in Fig. 5b by applying fracture mechanics. 

As Fig. 4 shows, at point P the curvature is not large, but torsion 
reaches a peak because the cilium is being twisted in the transverse 
plane. This is shown in the insets in Fig. 5a. Inset P shows that at 
point P the neutral plane of stress undergoes a reversal in com- 
pressive (C) and tensile (T) stresses. The two ends of the cilium in 
inset P also experience a reversal in the sign of the applied torque 
(t) between the cilium base and the distal point. Inset PI shows 
that, as a result, the two otherwise parallel circular cross-sections 
shear and open ajar due to the switching of the regions labeled C 
and T. 

The literature suggests that in the 9 + 2 axoneme the central 
microtubule pair is responsible for hardness control. It can also be 
surmised that motor protein molecules climb up the microtubule, 
find attachment spots, and lock 29 . This process is modeled in Fig. 5b 
and inset P2 in Fig. 5a. 
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Figure 6 | Collection of raw data and stages of processing of the temporal position of a cilium. (a) Distribution of the digitized cilium axial 
distance (um) from the base with beat phase (given by time step number); cilium of length 17 um. (b) Length distribution after smoothing along cilium 
length, (c) Length distribution after parametric fit of the closed orbits to smoothed axial data, (d) Side view of the cilium in the plane of the power stroke 
during the beat cycle from the processed data, (e) Top view of the cilium beat cycle from the processed data, (f) Cilium positions in the propulsion 
(x-y) plane (the red dot shows the base location; the base cilium is shifted horizontally to avoid clutter). 



The circled part of inset P2 shows the central microtubule pair 
{CM al4 and CM i, 14 ) and the CBL (CB P14 ). The presence of many 
such CBLs along the central microtubule pairs increases the second 
moment of area J ss , increasing the property G } ss , which is respons- 
ible for resistance to torsion, where G is the modulus of torsional 
rigidity of the material. (The CBLs are similar to crosslinks in 
polymers where they are closer when unstressed and pulled apart 
when stress is applied, thus straightening the polymer chain. The 
model of CBL attachment/detachment at P should apply to each of 
the microtubule pairs, and not just to the central pair of the 9 + 2 
axoneme.) 

In our model, when torsion at P reaches the maximum value, the 
CBL attachment cracks; this is shown in the boxed part of inset P2. 
The CBL at location P at TS14 is the proverbial "last straw" — when 
the CBL breaks, there is no hardness left for the power stroke to 
continue. Inset P2 shows how a crack develops (the sphere is the foot 
of the CBL) at the sites where the CBL is attached. The CBL motor 
protein and the host microtubule site are of dissimilar materials 
conformationally held in place. 

One expects cracks to develop where the CBL is attached to the 
microtubule and not at its apex because the microtubule attachment 
site has to be vacated fully for the new beat cycle to resume accepting 
a new CBL attachment. Our modeling considers the CBL cracking at 
the attachment site at TS14 at point P. The model includes CBLs all 
along the length of the microtubules, but only the one at point P 
cracks at TS14 at the contact site and needs replacing. 

One end of the protein is known to attach to one microtubule 
and the other end to another microtubule, creating a CBL structure 30 . 
An electrostatic model of the CBL attachment and detachment 
(although not in the context of the Paramecium's cilium beating) 
has been developed 29 . Here, fracture mechanics are applied to model 
the detachment. The CBL-microtubule contact is modeled to be 
brittle (like glass), rather than plastic (like metal); in other words, 
Griffith's law 31 applies and Irwin's 32 does not. Griffith's law states that 
the product of fracture stress and the square root of the crack length 
remain constant as the crack progresses until the difference between 
the surface and elastic energy reaches the peak value, at which point 
failure occurs. 



C 

Griffith's relationship may be applied as follows: oy = — =, 

l2Ey , ' V« 

C= \ , where oy is the stress at fracture, a is the crack length, 

V n 

C is a constant, E is Young's modulus of elasticity (= 2 GPa for 
polyethylene, = (10 to 30) GPa for myosin 33 ), and y is the surface 
energy density (= 0.72 J/m 2 for water, 1.0 J/m 2 for glass, and s 2 J/ 
m 2 for polymers; bulk cilium is modeled as a water-like polymer). 
Other parameters are as follows: equilibrium spacing of the CBL = 
10 nm, and the radius of the actin filament = 5.5 nm 29 . If we assume 
the value of a to be (5.57i) nm (hemispherical arc length), and E ( = 
0.50 GPa) of human tendon 34 to be applicable (A, the hemispherical 
area 27i(5.5 nm) 2 ), then the fracture stress (ay) would be 

= / 2(0 5x 10 )0.72 = x force required 

; V 7t 2 5.5xl0- 9 H 
to fracture one of the foot anchors of the CBL = 2.6 X 10~ 8 N = 
26 nN = 26 X 10 3 pN. 

Earth's gravity exerts a force of 10 - 100 pN on cells in bio- 
physical systems. In the artificial gravity experiments with parame- 
cia, the animals respond to gravity and regulate their swim speed 
primarily due to the buoyancy of the cell 35 . In a thermal cum elec- 
trostatic modeling of the motion of a single motor protein (myosin) 
molecule, a maximum CBL force of 5 pN has been assumed to 
apply 29 . So, a CBL detachment force of 26 X 10 3 pN is too high, 
and 5 pN is more reasonable. 

On the other hand, Cordova et al 29 state that (1) the attachment 
and detachment sites of the CBL differ by the presence of ATP when 
detached; and (2) the "binding of ATP weakens the attachment of the 
CBL to the fiber." In linear elastic fracture mechanics, Griffith's 
theory states that free energy is given by the difference between the 
surface energy and the elastic energy, and it increases with the crack 
length. Failure occurs when the free energy attains a maximum value 
at a critical crack length beyond which the free energy decreases by 
increasing the crack length. The free energy in Griffith's theory may 
be corrected due to the hydrolysis of ATP at the detachment site 
when the chemical energy is released. One can back- calculate what 
this should be as follows. The force (F) that leads the CBL to detach = 
[A, the hemispherical attachment area X oy] = 27t(5.5 X 10~ 9 ) 2 oy(N) 
= 190 X 10-' 8 oy (N). For this force (F) to be 5 pN, oy = 5 X 10~ 12 / 



SCIENTIFIC REPORTS | 3 : 1956 | DOI: 1 0. 1 038/srepOl 956 



7 



(1.876 X 1(T 15 ) = 2.665 X 10 3 N/m 2 = 2.665 kPa. Assuming, as 
before, £ = 0.5 X 10 9 N/m 2 , the surface energy density (y) would 
have to be (1 J = 1 Nm) 3.856 X 1(T 10 J/m 2 = 38.56 nj/m 2 . 

The ATP hydrolysis has to convert the chemical energy at the 
detachment site to make the surface energy density = 38.56 nj/m 2 
for the force £ required (produced practically at a point) to detach the 
CBL to be 5 pN. The total energy ( = surface energy - strain energy = 
Gf 2 ntx 2 L 



(2£) 



= yai, where L is the third dimension of the crack, 



and strain energy is the elastic energy released in the fracture area due 
to the fracture; the strain energy comes from the torsion) spent per 
detachment is 75 X 1(T 24 J = 75 yj. 

Attributing the total energy to the chemical energy, the chemical 
energy spent is estimated as follows. One mol has 6.02 X 10 23 mole- 
cules; the energy released per ATP mol is 45.6 kj. Therefore, the 
energy released per molecule of ATP is 76 X 10 3 yj, which can power 
1000 CBL detachments. For 3000 cilia in one Paramecium, three 
ATP molecules are needed per beat cycle to detach one CBL in each 
cilium; so, at 14 Hz, 42 ATP molecules are needed per second (fewer 
are needed if the absolute viscosity drops). 

A full modeling of the ATP hydrolysis process is needed to quan- 
titatively complete the evaluation of the mechanism. The hardness 
control model should include the mechanisms of softening at the 
detachment point due to the presence of ATP, subsequent enhanced 
bonding due to ATP hydrolysis, the effect of torsion, and sliding 
between the microtubules. 

Discussion 

Why does torsion reach its peak value in the middle of the cilium? If 
torsion reaches its peak value (and the CBL breaks) too near to the 
base, then too much energy (and time) would be required to 
straighten the cilium; on the other hand, if torsion reaches its peak 
value (and the CBL breaks) too far distal, then drag would be higher 
during the return stroke. The middle location is therefore a com- 
promise for locating the CBL, the connecting spring. Until the CBL is 
mended, the cilium would act as a two-link pendulum, which is 
known to be chaotic. But instead, an initial-condition-independent 
locked-in trajectory can ensure lowest energy consumption 5 . (This 
trajectory would be different in the cilia of different animals.) To 
achieve that, a Strouhal number defined as the ratio of the product 
of beat frequency and arc length of excursion divided by the forward 
speed of the Paramecium would have to be held to very high pre- 
cision (about eight decimal places). The agreement in Fig. 2 is 
remarkably good, indicating that this indeed may be the case. (The 
cilium beat is a resonant oscillation with little friction and amplified 
oscillations — see eq. (1)). To ensure precisely that the mass ratio of 
the cilium on each side is maintained, the broken CBL (or its recept- 
acle) may be of somewhat specialized construction. We also specu- 
late that this CBL acts as a multiplexer and demultiplexer 
(mux-demux); since the available time, the number of pathways 
for information, and the bandwidth are limited, to maximize 
throughput along the cilium, it is optimal to locate the switches near 
the middle. The switch idea is not required by the present model, but 
it serves to kill two birds with one stone. In the multi-systemic con- 
text, the modeled autonomic motion oscillator is coupled to numer- 
ous other oscillators in the animal. 

If we add a damping term to the simple oscillator equation given 
earlier, the expression for the unforced vibration of a damped 
pendulum is JO + 2 c T 0 + z k0 = 0, where 0 is the torsional angle, / 
is the moment of inertia of the pendulum, and c T , T k are torsional 
damping and spring coefficients 36 , respectively. The moment of in- 
ertia 37 of a rod-like pendulum of length L about its base is / = 
<nd 2 L\ 

mL 2 P l 



This number was obtained using values of p = 1000kg/m 3 , 
d = 0.25 um, and L=17 urn. The torsional spring coefficient ( z k) 

Nm 

for a kinocilium in a turtle's utricle 37 is 1.77 + 0.47 * 10~ 16 — -. 

rad 

The undamped natural frequency 38 of this torsional oscillator is 

Fk 

co„ = < / — = 46900 r jf c = 7.47 kHz, which is broadly in line with the 



figure of 16 kHz characteristic of the RLC model described below, 
but still far in excess of the 14 Hz found experimentally. The damped 



natural frequency of an oscillator 38 is coj = ft 
damping ratio (() is given by the expression £ = 



The damping 



- = — pd 2 ^ = 8.04 * W 2b kem 2 , where p and 
3 3 12 

d are the volumetric density and diameter of the rod, respectively. 



ratio required to yield a damped natural frequency of 14 Hz for a 
torsional oscillator consisting of the cilium is £~0. 999998, a value 
vanishingly close to the critical damping ratio of unity. 

Considering ionic circuits, to straighten the cilium during the 
return stroke, the Paramecium fuel cell provides energy that is used 
during the power stroke. The energy in a given amount of ATP is 
used to oscillate the cilium at the frequency it does at a given tem- 
perature. The presence of the ATP turns the cilium into an RLC 
circuit and the oscillatory energy transfer between the inductive 
(L) and the capacitive (C) components occurs optimally at the res- 
onant frequency. The load, which varies with temperature, is mod- 
eled by resistance _R, the current reaches the maximum value at the 
resonant condition, and the ATP consumption can be estimated. If C 
is 10 nF and L is 10 nH, the natural resonance frequency (/„ = (1/ 
(2ti)) * [1/^(LQ]) is 16 MHz. This is X 10" of 14 Hz, the beat 
frequency of the cilium in water at 20°C Spectral peaks at 
370 kHz (attributed to the myosin thick filament) and a relaxation 
time of 5 us (attributed to the CBL motion) have been observed in a 
cilium 3 '. If more CBLs link the microtubule pair, increasing C and L 
to 10 uF and 10 uH, respectively, the resonant frequency would be 
16 kHz, which is X 10 3 of 14 Hz — a value (roughly) closer to that 
from the structural consideration and the spectral peak of 370 kHz. 
Parametric oscillators with second-order nonlinear interactions can 
produce lower frequencies from two higher frequencies by differen- 
cing their frequencies; however, there is no evidence that this is 
happening. An explanation in the context of self-regulation of how 
C is held to high precision in order to yield a low value of damped 
frequency is more amenable to the olivo-cerebellar model given 
earlier. 

In the context of messaging mechanisms 7,40 , ATP oscillators cou- 
ple with the Ca 2+ -oscillator to regulate the CBL architecture, which 
alters the beat frequency. The precision of hardness can be controlled 
one CBL link at a time — a precision beyond the reach of current 
engineering. The precision of the damping ratio control is similar 
to the precision of the Strouhal number control of oscillation in the 
vortex-based flapping fin propulsion of swimming and flying ani- 
mals 5 , whose Reynolds number is higher. The, the self-regulation 
mechanism of control of the damping ratio to yield resonant amp- 
lification of actuator excursion possibly remains invariant with 
Reynolds number. 

Consider the efficiency of a Paramecium. The input power to the 
entire Paramecium for all its activities may be estimated as follows. 
Measurements of the oxygen respiration rate of a colony of 
Paramecium caudatum* 1 indicate a value of 6 (1IO2/IO 4 organism 
hr. At standard temperature and pressure, this is equivalent to 
o r = 2.36 * 10" 15 mol/org. sec. The energy yield from ATP due to 
aerobic respiration 42 in mitochondria is £ = 2880 kj per 6 mol O2, 
meaning that each Paramecium consumes approximately Wo 2 = 
o r £=1.13nW. 

The output thrust power may be estimated as follows. Since Re is 
very small, a Stokes flow approximation may be used to estimate the 
power needed by a Paramecium to move. The drag force 

Fd = — CopU 2 A, where Cd,P,U, andA are the drag coefficient, fluid 
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density, forward speed, and frontal area, respectively. We will assume 
the body length (B L ) of the Paramecium to be =200 urn, and the 
forward speed to be fc Bi/ sec, with its shape being approximately a 

2 : 1 prolate spheroid such that A = n (^^j j 4- The drag coefficient 

(C_>) can be assumed to vary as the reciprocal of Re, such that 

Co~C* /Regi. The drag force is Fj}= — C pv k Bi . The power 

consumed by this drag (output power) is Wd^FqU, such that 

Wp = —C* pv ic 2 Bi 3 =2.7 pW, where the above value assumes 

£ L = 200 um, k=12, p = 1000 kg m _3 ,v=l. * 10~ s m 2 s" 1 , and 
C*=24 (sphere). Taking our earlier estimate of cilium thrust effi- 
ciency of 0.125, and considering the Paramecium's numerous 
appended cilia, we lower the thrust efficiency of the organism to, 
say, 0.084 43 . The power input to thrust from oxygen can be obtained 

from i-] D = =-^- =0.084, with W T(h = 0.032 nW. A very small 

Wro 2 

amount (3%) of the total oxygen power is used for thrust suggesting 
excellent cost control. 

Note that the sudden wilting of the viscosity-dominated cilium 
motion, triggered by CBL breakup, is effectively a symmetry breaking 
that lowers drag in the return stroke, yielding a net thrust. Here, 
thrust is produced by laminar streams". On the other hand, in in- 
ertia-dominated swimming animals, thrust is produced by reverse 
Karman vortex jets. Here, the hydrodynamic (thrust) efficiency is 
higher: a single penguin-like flapping fin has values of about 0.60 6 , 
and a sunfish operating two pectoral flexible flapping fins has a value 
of 0.40 (due to Lauder & Madden — see Bandyopadhyay and 
Leinhos 43 ); (this value has been reproduced in a robotic model whose 
propulsion power density matches that of a shark 43 ). If improvement 
in efficiency was a motivation for the evolution of larger swimming 
animals, then the symmetry breaking 44 of the drag vortex wake to a 
jet vortex wake was a significant milestone that accompanied the 
increase in the portfolio of motions. This made fin-wake coupling 
and acoustic localization feasible, leading to predator-prey interac- 
tions 5 . The critical Reynolds number of symmetry breaking of wake 
vortex may have coincided with a significant increase in the number 
of inferior- olive neurons in the swimming animal. 

In hydrodynamic models incorporating the detailed 9 + 2 axo- 
neme internal elastic structure, the cilium beating appears as an 
emergent property of the coupling of the fluid and the actuator 
architecture 45,46 . Similarly, in complete swimming and flying animals 
that propel based on vortex dynamics, the hydrodynamics, control, 
and sensing follow similar nonlinear self-regulating equations — the 
fin and wake flow couple, and the fin architecture emerges as elastic, 
absorbing energy from the wake at a certain phase and releasing that 
energy at some other appropriate phase 5 . The benefit of multi- 
systemic self- regulation is that if sensors, actuators, and their control 
simultaneously follow the same dynamic systems principles to 
remain in persistent synchrony with the environment, then they 
home faster on moving targets if they are assumed to have a 10- 
15% stereoscopic bias (handedness) in sensing and actuation 5 . 

The present model of building fontanels in otherwise hardened 
structures for hardness control offers an inspiration to the nonlinear 
engineering design of actuators operating near neutral equilibrium at 
an unprecedented level of precision. In a reversal of the role of the 
Ca 2+ -oscillator, the cilium can be similarly modeled as a nonlinear 
sensor. The attractor basin of the cilium's LCO in Fig. 2a is a vast 
accurate library of phase. The trajectory of the temporal states of the 
sensor in the basin, and of their speed of return to the LCO, can be 
used to determine the amplitude and phase of any external perturba- 
tion that disturbs the cilium receptor off the LCO. 

Basic ciliary biology is translating to biomedicine where a large 
amount of information is being gathered that needs to be under- 
stood 7 . The present work suggests a relationship of ciliary biology 



with olivo-cerebellar neuroscience, which has a vast set of bifurcation 
properties 47 . It may be that ciliary mechanisms can be modeled as 
coupled nonlinear oscillators, which are deterministic but unpredict- 
able due to dependence on initial condition. Perhaps there are weak 
disparate oscillators that control the main oscillator, and a loss of 
coupling is manifested as a disease. 

Methods 

Acquisition and processing of the cilium position data. The University of British 
Columbia (UBC) Biology Department has produced an animation of the cilium 
motion using reported high-speed stereo-microscopy pictures 4 . Therein, the cilium 
position, velocity, and acceleration do not visually exhibit any spatio-temporal kink 
(smooth gradients of time and length). These cilium positions were digitized, 
parametrically processed, and smoothed (Fig. 6). The cilium motion is given by the 
position vector (in pixels) of a space curve r(s,t) = {x(s,t),y(s,t), z(s,t)} (Fig. SI-2). The 
cilium is divided into 20 segments (initially of non-uniform spacing). We define 5 as 
the axial length from the base to the distal point (L) of the particular cilium segment; 
also, we define an orthogonal coordinate system (x, y, z) that represents the propulsive 
(or stream-wise), vertical, and spanwise directions, respectively, and t = time. Due to 
cilium bending, initially, the positions (in pixels) along the length were not equally 
spaced in the longitudinal and top-view planes. For this reason, the digitization had to 
be processed in the following manner. With cilium bending, the (y(s,t)) distances 
could be determined more accurately from the (x(t),y(t)) view (side view), and the 
(x(s,t), z(s,t)) position values could be determined from the (x(t), z(t)) view (top view). 
The initial length was obtained from the three coordinates obtained this way (Fig. 6a). 
From the polynomial fits at each time step, the cilium was digitized at equal length 
intervals (r(s, f)). The data were further smoothed along the orbits (Fig. 6b). 
Subsequently, a three-dimensional parametric fit (Fig. 6c) was applied 28 . The result 
was a smoothed position distribution both along the cilium lengths at each time step 
and also of their rotational orbits obtained from position vectors at equispaced 5 
locations. The improvement in accuracy due to axial and rotational track smoothing 
of the cilium position vector data is given in tables SI-1 through SI-3. The cilium 
length of 1 7 um is maintained at all time steps; the local lengths are also maintained at 
all positions along the length of the cilium at all time steps. 

Source of unprocessed experimental position data of cilium. Fig. lb: measurements 
are sourced from Ref. 4. 

Figs. 2-6: original cilium position data as sourced and adapted from Ref. 4 by 
Prof. J. Berger of the University of British Columbia Biology Department; 
http://www.zoology.ubc.ca/courses/bio332/flagellar_motion.htm. 

Added Mass. Due to fluid-structure interaction, in a long mechanical cilium of fixed 
hardness and square section, added mass effects may account for a part of its low 
values of resonant frequency 48 . But, these effects are less in a smaller cilium of circular 
section which may imply some kind of hardness control in the biological cilium. An 
explanation in the context of self- regulation of how £ is held to high precision in order 
to yield a low value of damped frequency is more amenable to the olivo-cerebellar 
model given earlier. 
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